ENTAT10N  PAGE 


Form  AofitVY^a 

am  No.  0704-0703 


«at*a  to  ntrtqe  ‘  'ouf  oar  'ncatte.  ncuaing  tn«  n m*  iof  rw»?snnq  namcjotn.  iMfOttn*  erteivg  uu  vsurrsi. 
I  Tr**>rq  :n«  ton««ion  at  "(onnation  '♦no  com*«*o  roasrott'fl  aw  owosa  fsoowo  or  *»,  otrser  nsstt  a*  ,f'i 
t oixoan  to  «*  w i*qtoo  tacoouarton sanicw. Oirottcrata Aa  iniormaoon  Ocaraojfo <na aacara.  J.*. trtorwrt 
r  Ottico  o'  mo  Suaoat.  **awoan  ProrKt  ;37TSA-«itt).  /reywxjxa*.  OC  .asai 


CRT  DATE 


4.  TITLE  AMD  SU8TTOE 

STUDIES  IN  STATISTICAL  SIGNAL  PROCESSING 


4.  AUTHORS) 

Thmas  Kailath 


.  PERlKJaRSSftS  ORGANIZATION  WAME(S)  AND  ADORSSS(ES) 

Information  Systems  Laboratory 
Department  of  Electrical  Engineering 
Stanford  University 
Stanford,  CA  94305-4055 


3.  REPORT  TYPE  AND  OATES  COVERED 

Final  Report,  1  Jul  88  to  30  Jun  90 

S.  FUNDING  NUM3SR3 

AFOSR-88-0327 
61102F  2304/^6 


AF0SR.TR. 


6.  PERPORMIH6  ORGAMIXATIQM 
REPORT  NUMBER 


<>"  Of)  5 5 


ITsFONSOIUNG /  MONITORING  AGENCY  NAME(S)  AND  ADOSESS<IS) 

AF0SF./M 

SSrtMM  20332-6448 


10.  SPONSORING  /  MONITORING 
AGIKC?  SSPOST  NUM3ER 


I  iC 


AFOSR-88-0327 


12a.  OttTmaunON/  AVA&ABJUTY  STATEMENT 

y.  t  ip  i 

for  DUDi  aC  >  - 

iPM%r4«l«UDll“lt°4' 


13.  ABSTRACT  (Mixtmum  200  wotxH) 


This  is  a  final  report  on  grant  AFOSR-88-0327,  for  the  period  of  July 
1,  1988  to  June  30,  1990.  Section  1  provides  a  brief  overview  of  our 
work,  while  Section  2-5  describe  in  some  detail  our  recent  results  on 
efficient  factorization  of  structured  matrices,  recursive  layer 
peeling,  recursive  state-space  synthesis  and  wavelet  representations. 


14.  SUBJECT  TERMS 


IB.  PRJCT  CODS 


r?T  sEajwTTOMStfScJrioiri 

OP  REPORT 

UNCLASSIVIED 


18.  SECURITY  CLASSIFICATION 
OP  THIS  PAGE 

UNCLASSIFIED 


9.  SECURITY  OASSlftCAnCN 
OP  ABSTRACT 

UNCLASSIFIED 


|  20.  UASTATION  OF  ASSTOACT 


SAR 

Standard  Form  298  (Rev  2  39) 

■"♦ACnwe  Dv  AAiV  •  39-  1  © 


Department  of  The  Air  Force 


FINAL  REPORT 

for 

July  1,  1988  to  June  30,  1990 

Contract  AFOSR  88-0327 
Task  No.  2304/ A6 


Accession  For 


NTIS  GRA&I 
DTIC  TAB 
Unannounced 
Justification. 


□ 


By- 


Distribution/ 


Availability  Codes 


Dlst 


to 


Avail  and/or 
Special 


Studies  In  Statistical  Signal  Processing 


Thomas  Kailath 
Information  Systems  Laboratory 
Department  of  Electrical  Engineering 
Stanford  University 
Stanford,  CA  94305-4055 


The  views  and  conclusions  contained  in  this  doc¬ 
ument  are  those  of  the  author  and  should  not  be 
interpreted  as  necessarily  representing  the  official 
policies  or  endorsements,  either  expressed  or  im¬ 
plied,  of  the  Air  Force  Office  of  Scientific  Research 
of  the  U.S.  Government. 


STUDIES  IN  STATISTICAL  SIGNAL  PROCESSING 


This  is  a  final  report  on  Contract  AFOSR-88-0327,  for  the  period  of  July  1,  1988 
to  June  30,  1990.  Section  1  provides  a  brief  overview  of  our  work,  while  Sections  2-5 
describe  in  some  detail  our  recent  results  on  efficient  factorization  of  structured  matrices, 
recursive  layer  peeling,  recursive  state-space  synthesis  and  wavelet  representations. 

1.  INTRODUCTION 

The  primary  objective  of  our  research  is  to  develop  efficient  and  numerically  stable  al¬ 
gorithms  for  nonstationary  signal  processing  problems  by  understanding  and  exploiting 
special  structures,  both  deterministic  and  stochastic,  in  the  problems.  We  also  strive 
to  establish  and  broaden  links  with  related  disciplines,  such  as  cascade  filter  synthesis, 
scattering  theory,  numerical  linear  algebra,  and  mathematical  operator  theory  for  the 
purpose  of  cross  fertilization  of  ideas  and  techniques.  These  explorations  have  led  to  new 
results  both  in  estimation  theory  and  in  these  other  fields,  e.g.,  to  new  algorithms  for 
triangular  and  QR  factorization  of  structured  matrices,  new  techniques  for  root  location 
and  stability  testing,  new  realizations  for  multiple-input/multiple-output  (MIMO)  trans¬ 
fer  functions,  and  new  recursions  for  orthogonal  polynomials  on  the  unit  circle  and  the 
real  line  as  well  as  on  other  curves. 

For  several  years,  the  guiding  principle  in  these  studies  has  been  the  concept  of  gen¬ 
eralized  displacement  structure  (Lev-Ari  and  Kailath  (1986)),  which  generalized  and 
subsumed  our  earlier  work  on  Toeplitz-oriented  displacement  structure  (Kailath,  Kung 
and  Morf  (1979);  see  also  Lev-Ari  and  Kailath  (1984)).  A  related  notion  of  displace¬ 
ment  structure  has  also  emerged  from  recent  work  of  Heinig  and  Rost  (1984,1987)  of 
East  Germany.  While  they  are  aware  of  our  work,  and  make  some  attempts  to  relate  to 
it,  their  approach  and  methodology  are  significantly  different  from  ours.  In  particular, 
they  focus  only  on  the  problem  of  inversion  of  structured  matrices  via  algebraic  methods, 


while  our  work  has  primarily  addressed  triangular  factorization  of  such  matrices,  and  our 
approach  is  based  on  a  generating  function  characterization  of  matrices.  The  triangular 
factorization  problem  is  in  many  senses  more  fundamental  than  inversion,  and  has  more 
consequences  for  signal  processing,  linear  algebra,  operator  theory  and  other  fields.  In 
fact,  recently  we  were  able  to  reduce  the  inversion  problem  for  structured  matrices  to 
the  factorization  of  certain  block-matrices  with  structured  blocks  (see  Chun  and  Kailath 
(1989)).  This  result  also  confirms  and  clarifies  an  earlier  observation  [see  Lev-Ari  and 
Kailath  (1984)]  on  the  relation  between  efficient  inversion  and  efficient  factorization  of 
structured  matrices:  only  some  of  the  structured  matrices  that  admit  an  efficient  factor¬ 
ization  procedure  can  also  be  efficiently  inverted. 

The  generating  function  approach  also  suggests  a  natural  system-theoretic  interpre¬ 
tation  of  the  theory,  which  allows  a  study  of  various  problems  in  system  theory,  such  as 
minimal  realization,  Pade  approximation,  control  design,  and  a  variety  of  root  distribu¬ 
tion  (stability)  problems  for  polynomials.  Perhaps  the  most  prominent  system-theoretic 
aspect  of  our  efficient  factorization  techniques  is  that  they  can  be  interpreted  as  recursive 
identification  procedures  for  certain  lossless  cascade  models.  For  instance,  the  classical 
Schur  algorithm  is  also  a  procedure  for  step  by  step  identification  and  ‘peeling’  of  the 
layers  of  a  transmission-line  with  a  piecewise  constant  characteristic  impedance  [Bruck- 
stein  and  Kailath  (1987)].  Such  layered  models  have  been  used  for  quite  some  time  in  oil 
exploration  and  in  marine  seismography.  They  involve  two  scalar  signals  propagating  in 
opposite  directions;  consequently,  the  characteristics  of  the  model  can  be  captured  by  a 
single  scalar  input-scattering  function.  Schur’s  original  formulation  of  his  algorithm  was, 
in  fact,  in  terms  of  this  scattering  function. 

We  have  recently  begun  to  extend  our  methods  to  multichannel  cascade  models ,  which 
involve  involve  multiple  signals  propagating  both  in  the  forward  and  in  the  opposite 
(backward)  direction.  Since  such  models  are  represented  by  matrix  scattering  functions , 
it  would  seem  that  the  corresponding  layer-peeling  procedures  need  to  be  rederived  in 
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matrix  (or  block)  form.  This  is  certainly  possible  (see,  e.g.,  Delsarte,  Genin  and  Kamp 
(1979)),  but  results  in  the  introduction  of  computation-intensive  matrix  operations,  such 
as  the  evaluation  of  matrix  square  roots.  In  contrast  to  this  approach,  we  have  succeeded 
in  obtaining  layer-peeling  procedures  that  involve  only  elementary  (2  x  2)  circular  and  hy¬ 
perbolic  rotations,  and  therefore  require  only  scalar  computations.  We  have  achieved  this 
by  incorporating  the  notion  of  modular  decomposition ,  which  has  been  used  in  the  past  to 
decompose  the  multichannel  Levinson  algorithm  (see  Sakai  (1982))  and  the  deterministic 
recursive-least-squares  lattice  filter  (Lev-Ari  (1987)).  Our  modularly-decomposed  filters 
can  be  implemented  in  pipelined  parallel  processing  hardware  (such  as  systolic  arrays), 
with  the  throughput  being  independent  of  the  number  of  channels  ( i.e .,  the  number  of 
forward  and  backward  signals  flowing  through  the  model). 

Moreover,  our  modular  formulation  makes  it  possible  to  extend  the  Schur  algorithm  to 
matrix  functions  with  poles  within  the  unit  disc,  which  arise  in  various  control  problems 
and  in  particular  in  model-order  reduction  with  Hankel  norm  (see,  e.g.,  Genin  and 
Ivung  (1981)).  Such  an  extension  would  have  been  impossible  in  the  block-formulation 
of  Delsarte,  Genin  and  Kamp  (1979);  and  it  would  be  quite  difficult  to  obtain  even  in 
the  ‘tangential’  formulation  of  Dewilde  and  Dym  (1981).  Moreover,  we  have  applied  the 
same  approach  to  extend  in  a  similar  way  also  the  Nevanlinna  algorithm,  which  includes 
Schur’s  algorithm  as  a  particular  case  (z'.e.  ,  all  extraction  points  are  at  the  origin). 

While  the  Schur  and  Nevanilnna  algorithms  synthesize  lossless  cascade  models  from 
their  input  scattering  function,  there  has  been  also  much  work  done  on  the  synthesis  of 
lossless  arrays  from  state-space  representations  (see,  e.g.,  Genin  et  al.  (1983),  Roberts 
and  Mullis  (1987)).  This  approach  involves  more  computation,  but  results  in  a  greater 
flexibility  in  choosing  the  array  configuration  (which  need  not  be  a  cascade).  In  addition, 
it  provides  convenient  analytical  tools  for  the  analysis  of  finite  precision  effects  such  as 
overflow  oscillations  and  roundoff  noise. 

We  have  recently  formulated  a  new  (recursive)  approach  to  lossless  cascade  synthesis 
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from  state-space  representations.  Our  approach  not  only  requires  significantly  less  com¬ 
putations  than  the  traditional  method  (which  is  based  on  conversion  to  balanced  form), 
but  also  leads  in  a  natural  way  to  a  matrix-domain  formulation  of  the  generalized  Schur 
algorithm  for  triangular  factorization  of  structured  matrices.  This  new  formulation  sub¬ 
sumes  most  previously  published  procedures  for  factorization  of  structured  matrices  and. 
in  particular,  the  matrix-domain  formulation  of  Chun  and  Kailath  (1989,  1990)  and  the 
transform-domain  formulation  of  Lev-Ari  and  Kailath  (1984,  1986).  Finally,  the  state 
space  approach  is  also  helpful  in  describing  the  relation  between  the  cascade  models  as¬ 
sociated  with  structured  matrices  (as  described  in  Lev-Ari  and  Kailath  (1986))  and  the 
models  associated  with  the  inverses  of  these  matrices.  These  ‘inverse’  models  are  the 
starting  point  for  the  derivation  of  the  Levinson  algorithm  and  the  Gohberg-Semencul 
formula  for  structured  matrices  other  than  Toeplitz. 

Another  application  of  the  state-space  approach  has  been  presented  by  Glover  (1984) 
in  the  context  of  model-order  reduction.  We  have  recently  begun  to  explore  the  possibility 
of  his  method  with  our  results  on  multichannel  Schur  algorithms  in  order  to  obtain  a 
computationally-efficient  procedure  for  determining  a  modular  realization  of  the  reduced 
order  filter. 

So  far  we  have  considered  only  (linear)  models  with  time-invariant  parameters.  In 
order  to  be  able  to  characterize  processes  with  stable  (i.e.,  bounded)  persistently  non¬ 
stationary  dynamics  we  need  to  allow  models  with  time-varying  parameters.  For  instance, 
processes  with  periodically-varying  correlation  (i.e.,  rt+piT+p  —  rt,s  for  all  t,s)  involve 
models  with  periodically-varying  coefficients.  One  way  to  overcome  such  complications  is 
to  consider  representations  cf  processes  as  linear  combinations  of  known  functions,  simi¬ 
lar  to  the  Karhunen-Loeve  representation.  However,  we  focus  on  basis  functions  that  are 
independent  of  the  statistics  of  the  process  in  consideration.  In  particular,  we  have  con¬ 
sidered  basis  functions  derived  from  the  theory  of  wavelets  (see,  e.g.,  Daubechies  (1988), 
Mallat  (1989)).  We  have  characterized  the  asymptotic  behavior  of  such  representations 
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with  increasing  level  of  resolution  and  have  established  (so  far  only  for  periodically- 
correlated  processes)  the  correlation  structure  of  the  representation  coefficients. 

2.  FACTORIZATION  OF  STRUCTURED  MATRICES 

Our  early  work  on  factorization  and  inversion  of  Toeplitz  and  close- to-Toeplitz  matrices 
lea  us  to  the  observation  that  for  certain  matrices  the  displacement  matrix 

VZR  :=  R  -  ZRZ*  ,  Z  =  [Silj+1\lj=0 

has  low  rank.  Notice  that  Z  has  unity  elements  on  the  first  subdiagonal  and  zeros 
elsewhere.  Consequently,  the  displacement  matrix  VZR  is  the  difference  between  the 
matrix  R  and  the  matrix  ZRZ*  obtained  by  displacing  R  one  step  along  the  main 
diagonal.  In  particular,  the  displacement  rank  (i.e.,  the  rank  of  VZR)  is  2  for  both 
Toeplitz  matrices  and  inverses  of  Toeplitz  matrices.  We  have  shown  in  previous  work 
(largely  supported  by  AFOSR)  that  the  displacement  concept  is  a  key  tool  for  developing 
fast  algorithms  of  many  kinds,  including  factorization  and  inversion  of  Toeplitz  and  near- 
Toeplitz  matrices,  as  well  as  fast  (generalized  Levinson  and  Schur)  algorithms  for  solving 
linear  systems  with  such  coefficient  matrices.  Not  surprisingly,  these  results  led  naturally 
to  cascade  orthogonal  structures  for  the  prediction  of  nonstationary  processes  (Lev-Ari 
and  Kailath  (1984)).  We  have  also  found  that  the  same  concept  is  tightly  connected 
to  the  more  general  problem  of  cascade  filter  synthesis  in  network  theory  and  digital 
filtering  as  well  as  to  a  variety  of  inverse  scattering  problems  (some  references  are  Rao 
and  Kailath,  (1984,  1985),  Bruckstein  and  Kailath  (1987),  Lev-Ari  (1988)). 

Later  we  extended  the  concept  of  displacement  structure  to  a  very  broad  family 
of  structured  matrices,  including  Hankel  matrices  and  their  inverses,  sums  of  Toeplitz 
and  Hankel  matrices  and  several  others  (Lev-Ari  and  Kailath,  (1986)).  The  generalized 
displacement  of  a  matrix  R,  is  defined  as  d( Z,  Z)R  where 

d( A,B)R:=  £  4,iAfcR(B*y  ,  (1) 

k,l=  0 
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and  the  asterisk  (*)  denotes  Hermitian  transpose  (complex  conjugate  for  scalars). 

The  concept  of  displacement  structure  and  its  properties  are  more  conveniently  de¬ 
scribed  in  terms  of  generating  functions.  The  generating  function  of  a  matrix  R  is  a 
power  series  in  two  complex  variables,  viz., 

R(z,w)-.=.[  1  2  z1  ...]R[1  w  w 2  ...]*  (2) 


The  displacement  d( Z,  Z)R  of  a  matrix  has  the  generating  function  d(z,  w)R(z,  w),  where 

N 

d(z,w)=Y,  dk,izk(w’)1  (3a) 

k,l= o 


Thus  the  generating  function  of  a  Hermitian  matrix  with  a  displacement  structure  has 
the  form 


R(z,  w)  = 


G(z)  J  G‘(w) 


(3  b) 


d(z,w) 

where  J  is  any  constant  nonsingular  Hermitian  matrix.  The  triple  {d(z,w),G(z),  J} 
is  called  a  generator  of  R(z,w),  since  it  uniquely  determines  the  generating  function 
R(z,w). 


We  have  extended  our  previous  work  (Lev-Ari  and  Kailath  (1984))  on  efficient  factor¬ 
ization  of  matrices  with  displacement  structure  to  accommodate  the  generalized  displace¬ 
ment  d( Z,Z)R,  and  we  have  shown  (Lev-Ari  and  Kailath  (1986),  Lev-Ari  (1989)),  that 
efficient  factorization  of  R  is  possible  if,  and  only  if,  there  exist  power  series  a(z),  fi(z) 
(with  arbitrary  radii  of  convergence)  such  that 


d(z,w)  =  a(z)ot*(w)  —  f3(z)/3*(w)  .  (4) 

To  obtain  the  factorization  of  R  one  has  to  propagate  the  recursion  (with  G0(z )  =  G(z)) 

(z  -  ( i)Gi+1(z )  =  G;(z)0,'(z)  i  =  0, 1, 2, . . .  (5) 

where  ©{(z)  is  specified  in  terms  of  G,((,)  (see,  e.g.,  Lev-Ari  and  Kailath  (1986),  Lev-Ari 
(1989)).  The  lossless  (pd-<?)- port  ©1(2)  can  be  decomposed  into  a  constant  (memoryless) 
part,  which  has  p  +  q  inputs/outputs  and  a  scalar  (single-input/single-output)  dynamic 
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part  (Lev-Ari  (1989)).  This  decomposition  is  particularly  simple  when  d(0,0)  ^  0,  as 
described  in  Fig.  1. 


Figure  1.  Decomposition  of  ©;(z)  (p  =  1  =  q,  £  =  0). 

The  standard  choice  of  the  extraction  points  {(,},  i.e.,  Q  =  0,  produces  triangular 
factorizations;  other  choices  can  be  useful  in  root-location  and  filter  synthesis  procedures 
(see,  e.g.,  Deprettere  and  Dewilde  (1980),  Vaidyanathan  and  Mitra  (1984)).  This  algo¬ 
rithm  requires  0(n 2)  computations  to  factor  a  structured  n  x  n  matrix  R  in  the  form 
R  =  LDL*,  in  contrast  to  the  conventional  LDL *  algorithm  which  requires  0(n3)  oper¬ 
ations  to  factor  an  arbitrary  n  x  n  matrix.  The  z-th  element  of  the  diagonal  matrix  D 
and  the  z-th  column  of  the  lower  triangular  matrix  L  are  determined  by  the  coefficients 
of  the  power  series  expansion  of  Gi(z). 

3.  RECURSIVE  LAYER  PEELING 

The  fundamental  factorization  procedure  (5)  for  structured  matrices,  viz., 

(z-C,-)G,-+1(2)  =  Gt-(2)0^) 

is,  at  the  same  time,  also  a  layer-peeling  procedure.  Starting  with  Go(z ),  which  we  can 
interpret  as  boundary  data  for  a  layered  medium,  we  identify  an  elementary  layer  with 
chain-scattering  matrix  0o(^)j  then  “peel”  it  off  to  obtain  G\(z),  the  boundary  data  for 
the  rest  of  the  medium  (with  the  first  layer  removed),  and  repeat  the  same  procedure 
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again  and  again.  Such  iayer-peeling  recursions  have  been  used  in  cascade  filter  synthesis 
(see,  e.g.,  Dewilde,  Vieira  and  Kailath,  (1978);  Vaidyanathan  and  Mitra  (1984)),  in 
inverse  scattering  (Bruckstein  and  Kailath  (1987)),  zero-location  (Lev-Ari,  Bistritz  and 
Kailath  (1987)),  and  model-order  reduction  (see,  e.g.,  Genin  and  Kung  (1981)). 

The  classical  work  of  Schur  (1917)  forms  the  basis  for  much  of  the  subsequent  work  on 
layer  peeling  procedures.  Schur’s  algorithm  associates  a  sequence  of  so-called  reflection 
coefficients,  all  with  magnitude  bounded  by  unity,  with  every  passive  scattering  function, 
i.e.,  a  function  f(z)  that  is  analytic  and  bounded  by  1  in  the  unit  disc.  In  particular, 
if  f(z)  is  an  all-pass  function,  which  means  that  |/(ej9)|  =  1  for  all  0 ,  then  Schur’s 
algorithm  produces  a  finite  sequence  of  reflection  coefficients  {k,  ;  0  <  i  <  n),  where 
\kn\  =  1  and  \k{\  <  1  for  0  <  i  <  n  —  1. 

Another  property  of  the  algorithm  is  that  starting  with  a  passive  scattering  function 
it  generates  a  sequence  of  such  functions.  This  is  the  essence  of  the  layer  peeling  method: 
a  single  step  of  the  Schur  algorithm  applied  to  a  passive  medium  leaves  a  medium  with 
the  same  property,  which  makes  it  possible  to  apply  the  same  step  again  and  again.  A 
single  step  of  the  Schur  algorithm  corresponds  to  the  removal  (or  peeling)  of  an  elementary 
lossless  two-port.  Thus,  the  algorithm  produces  a  discrete  transmission-line  model,  whose 
input  scattering  function  is  f(z)  (Fig.  2). 


i 


m 


•  • 


Figure  2.  Transmission-line  model  associated  with  the  Schur  algorithm. 


In  addition  to  the  recursive  characterization  of  passivity  via  the  constraint  on  the 
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magnitude  of  the  reflection  coefficients,  Schur  also  introduced  an  operator-norm  charac¬ 
terization  of  passivity:  he  proved  that  for  any  function  f(z)  that  is  analytic  in  the  unit 
disc  we  can  construct  an  infinite  lower-triangular  Toeplitz  matrix  whose  first  column 
consists  of  the  coefficients  of  the  power  series  expansion  of  /(z),  viz., 

ho  \ 

ft  fo  0 

L(/)  =  fo  f\  •  (6a) 


such  that 


sup  \f(z)\  =  ||L(/)||2  <  1 
N<i 


where  ||A||2  denotes  the  conventional  norm  of  a  matrix  A,  i.e., 

\\A\\,  :=  sup  (6c) 

*  IFII2 

and  |ja;||2  denotes  the  Euclidean  (£2)  norm  of  a  vector  x. 

Schur’s  algorithm  was  later  extended  by  Cohn  (1922)  to  functions  with  poles  in  the 
unit  disc ,  but  only  to  rational  “all-pass”  functions,  i.e.,  to  functions  f(z)  of  the  form 


f(z)  =  A 


p*(z) 


where  p#(z)  denotes  the  conjugate  reversal  operation,  viz., 

p*{z)  =  zd es  pW[p(l/2*)]*  (76) 

The  now  well-known  Schur-Cohn  test  associates  with  each  such  function  1  a  finite  se¬ 
quence  of  reflection  coefficients,  some  of  which  have  magnitudes  larger  than  1.  Moreover, 
it  has  been  shown  ( e.g .,  using  the  properties  of  Bezoutians  on  the  unit  disc)  that  the 
number  of  poles  of  f(z)  =  p#(z)/p(z)  inside  the  unit-disc  equals  the  number  of  singular 


Assuming  p(z)  has  no  zeros  at  z  =  0,  and  applying  the  algorithm  to  f(z)/\. 
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values  of  the  matrix  L(/)  that  are  larger  than  |A|  or,  equivalently,  the  number  of  negative 
eigenvalues  of  the  following  finite  rank  matrix, 

R:=|A|2I-L(/)L*(/)  .  (8) 

Rational  allpass  functions  of  a  given  degree  k  are  members  in  the  family  H™  of  all 
functions  with  k  poles  (or  less)  inside  the  unit  circle,  and  whose  magnitude  is  bounded 
on  the  unit  circle,  i.e., 

sup  [f(z)\  <  oo 

l*l=i 

It  turns  out  that  the  Schur-Cohn  algorithm  does  not  map  the  family  H into  itself,  ex¬ 
cept  when  k  =  0.  This  means  that  this  algorithm  does  not  admit  the  same  layer-peeling 
interpretation  as  the  standard  Schur  algorithm.  Nevertheless,  we  have  found  that  it  is 
possible  to  modify  the  Schur  algorithm  in  such  a  way  that  the  resulting  recursion  indeed 
maps  the  family  H into  itself  and,  therefore,  admits  the  same  layer-peeling  interpreta¬ 
tion  as  the  classical  Schur  algorithm.  Moreover,  the  layers  involve  only  elementary  (2  x  2) 
orthogonal  and  hyperbolic  rotations  (Ackner,  Lev-Ari  and  Kailath  (1990a)). 


P~ type  layer  (&,■  =  /,( o))  n-type  layer  (k{  =  i//*(o)) 

Figure  3.  Layers  for  the  modified  Schur  algorithm 

Furthermore,  our  modified  recursion  applies  to  every  function  f(z)  €  H £°,  and  not 
just  to  allpass  functions.  Each  layer  in  the  resulting  transmission-line  model  has  an  sign 
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or  ‘polarity’  (Fig.  3).  While  a  p-type  layer  (i.e.,  one  of  positive  polarity)  maps  into 
itself,  a  n-type  layer  maps  H%°  into  Hf£_ x  namely,  it  reduces  by  one  the  number  of  poles 
within  the  unit  disc.  Therefore,  the  number  of  n-type  layers  in  the  transmission-line 
model  that  is  generated  by  our  modified  algorithm  equals  the  number  of  poles  that  the 
function  f(z)  has  within  the  unit  disc.  This  is  the  key  idea  in  the  construction  of  efficient 
procedures  for  zero-location  and  in  the  solution  of  several  extension  problems  arising  in 
various  control  applications  (Ackner,  Lev-Ari  and  Kailath  (1990a)). 

In  the  layer-peeling  procedure  (5)  we  have  the  freedom  of  choosing  the  extraction 
points  {£,}.  The  flexibility  in  choosing  these  points  enables  us  to  optimize  the  number  of 
computations  in  the  algorithm  and  to  overcome  singularities.  The  standard  choice  £,  =  0 
corresponds  to  the  Schur  algorithm.  The  general  choice  corresponds  to  the  Nevanlinna 
algorithm  (Nevanlinna  (1929)).  In  (Ackner,  Lev-Ari  and  Kailath  (1990b))  we  extended 
our  results  from  (Ackner,  Lev-Ari  and  Kailath  (1990a))  to  the  Nevanlinna  algorithm 
and  showed  how  to  modify  the  recursions  in  the  case  of  meromorphic  functions.  The 
computational  procedure  corresponding  to  this  case  involves  Newton  Series  expansions , 
viz., 

03 

f(z)  =  £  fkVk(z) 

k=o 

where 

**(*)  ;=  n(*-o)  > 

3=0 

instead  of  the  MacLaurin  (i.e.  ,  power  series)  expansion  used  in  the  computational  form  of 
the  Schur  algorithm  (see,  e.g.  ,  Kailath  (1987)).  Applications  of  the  Nevanlinna  algorithm 
include  model-order  reduction  and  solution  to  several  other  interpolation  problems  arising 
in  control  and  signal  processing. 

Recursive  layer  peeling  for  a  medium  with  multiple  inputs  and  outputs  involves  a 
generalization  of  the  Schur  algorithm  to  matrix-valued  analytic  functions.  One  version 
of  the  matrix  Schur  algorithm  (Delsarte,  Genin  and  Kamp  (1979))  requires  hyperbolic 
matrix  rotations ,  which  are  computationally  expensive  since  they  involve  the  finding 
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of  the  square  root  of  (positive  definite)  matrices.  Moreover,  this  approach  cannot  be 
generalized  to  matrix  functions  with  elements  in  iJ|°,  because  now  it  would  involve 
square  roots  of  indefinite  matrices. 

An  alternative  approach  to  the  matrix  analytic  case  was  taken  by  Dewilde  and  Dym 
(1981).  They  extract  simpler  layers  than  in  the  method  of  Delsarte,  Genin  and  Kamp 
and,  as  a  resulf,  there  is  no  need  to  compute  square  roots  of  matrices.  This  computational 
procedure,  which  is  known  as  a  tangential  Schur  algorithm,  was  the  starting  point  for  our 
research  on  layer-peeling  methods  for  MIMO  (multiple-input/multiple-output)  systems. 


r  —  —  —  —  — 


Figure  4.  Single  layer  of  the  transmission-line  model  associated  with  the  multichannel 

Schur  algorithm  (p  =  3,  q  =  4) 
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We  have  shown  that  the  procedure  of  Dewilde  and  Dym  can  be  transformed  into 
an  equivalent  form  that  involves  only  elementary  (2  x  2)  rotations  (Ackner,  Lev-Ari  and 
Kailath  (1990c)).  In  fact,  for  p  x  q  matrix  scattering  functions,  the  peeling  of  each  layer  in 
our  version  of  the  algorithm  is  implemented  by  a  sequence  of  q  —  1  elementary  orthogonal 
rotations  {(jP;  0  <  k  <  p  —  1},  and  a  core  module ,  consisting  of  a  single-channel  layer 
(i.e.,  a  memoryless  lossless  two-port  with  reflection  coefficient  and  a  delay)  (Fig.  4). 
These  operations  can  be  easily  implemented  in  either  software  or  hardware. 

In  addition  to  requiring  significantly  fewer  computations  than  the  tangential  Schur 
algorithm,  our  version  also  serves  to  clarify  the  relationship  between  the  MIMO  case  and 
the  better-known  scalar  or  SISO  (single-input/single-output)  case.  The  MIMO  procedure 
differs  from  the  scalar  one  only  by  the  presence  of  elementary  orthogonal  rotations. 
Thus  both  procedures  share  the  same  core  module,  which  consists  of  a  single  elementary 
hyperbolic  rotation  and  a  (block-)  delay  element. 


Figure  5.  Structure  of  single  layer  in  the  non-analytic  case 
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Consequently,  we  are  able  to  show  that  our  version  of  the  tangential  Schur  algorithm 
can  be  modified  to  accommodate  scattering  functions  with  poles  within  the  unit  circle, 
and  that  this  modification  affects  only  the  core  module  (i.e.,  it  is  independent  of  the 
number  of  inputs  or  outputs).  This  means  that  each  layer  still  has  a  ‘sign’  or  polarity, 
as  in  the  SISO  case,  and  that  the  number  of  n-type  layers  coincides  with  the  number  of 
Smith-McMillan  poles  of  the  given  matrix  scattering  function.  Unlike  the  analytic  case, 
we  now  have  trivial  cells  (indicated  by  dashed  boxes)  both  in  the  upper  and  lower  parts 
of  the  diagram  (Fig.  5).  This  results  in  a  total  of  2P  possible  configurations  for  each 
layer,  in  contrast  to  a  single  configuration  in  the  (matrix)  analytic  case. 

A  single  layer  of  our  tangential  Schur  algorithm  involves  q  parameters:  one  in  the 
core  module  and  q  —  1  in  the  preceding  orthogonal  rotations.  In  comparison,  the  matrix 
Schur  algorithm  of  Delsarte,  Genin  and  Kamp  (1979),  viz., 

W  =  M,[I  -  F(z)F'( 0)]-1  [F(z)  -  F(0)]M,-> 

with 

m;m i  =  i  -  f(o)f*(o)  ,  m;m2  =  i  -  f*(o)f(o) 

which  corresponds  to  p  layers  of  the  tangential  algorithm,  has  the  same  number  of  pa¬ 
rameters  (p  x  q)  per  block-step. 

4.  RECURSIVE  STATE-SPACE  SYNTHESIS 

In  addition  to  the  input-output  approach  in  recursive  layer  peeling,  there  has  also  been 
a  considerable  interest  in  synthesis  of  lossless  models  from  state-space  descriptions.  For 
instance,  this  approach  has  been  used  extensively  in  orthogonal  synthesis  of  digital  fil¬ 
ters:  starting  with  a  given  state  space  model  for  an  allpass  filter,  one  first  obtains  an 
orthogonal  state  space  model  of  a  particular  (Hessenberg)  form  and  then  proceeds  to  ob¬ 
tain  a  realization  in  terms  of  elementary  (2  x  2)  lossless  cells  via  a  recursive  factorization 
procedure  (see,  e.g.,  Roberts  and  Mullis  (1987)).  Passive  transfer  functions  can  also  be 
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realized  in  this  fashion  by  embedding  them  into  allpass  filters  with  additional  inputs  and 
outputs.  This  can  be  accomplished  either  by  spectral  factorization,  as  in  the  method  of 
Rao  and  Kailath  (1984),  or  by  solving  certain  Ricatti  equations  (see,  e.g.,  Desai(1989)). 

To  be  more  specific,  denote  the  state  space  model  associated  with  a  given  lossless 
cascade  model  (such  as  the  one  described  in  the  previous  section)  by  {A0,Bo,  Co,  Do}. 
This  means  that  the  transfer  function  of  the  lossless  cascade  {0O(^),  ©1(2), . . . ,  0„(^)} 
can  be  expressed  in  the  form2 

0o(*)0i(*)  •  •  •  OnGO  =  Do  +  C0(z_1I  -  Ao)-1Bo  =  D0  +  zC0{I  -  *A0)_1B0  (9) 

Also,  it  can  be  shown  that  {Ao,  Bo,  Co,  Do}  is  J-orthogonal,  viz., 

A0  Bo  \  _  (  I  0 

Co  Do  J  ^  0  J 

where 

J  :=  diag{/p,  -Iq)  (106) 

Any  other  (minimal)  realization,  say  {A,B,C,D},  must  be^related  to  {A0,  B0,  C0,D0) 
via  a  similarity  transformation,  i.e., 

A  =  T~1A0T,  B  =  T-1B0  ,  C  =  C0T  ,  D  =  D0 

It  follows  that  the  J-orthogonality  relation  (10)  is  now  replaced  by 

a  bWr  oWa  b  \  _  /  R  0 
cdJ\ojJ\cdJ  \ 0  J 

where  R  :=  TT*.  In  particular  we  observe  that 

R  -  ARA*  =  B  JB”  (12) 

which  demonstrates  the  generalized  displacement  property  of  the  matrix  R.  This  ob¬ 
servation  has  led  Genin  et  al.  (1983)  to  propose  a  procedure  for  construction  of  lossless 

2 Recall  that  we  use  z  instead  of  z-1  to  denote  a  delay 


(  Ao  Bo  \  (  I  0  \  ( 

^C0  D0  /  \  0  J  J  ^ 
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cascade  model  by  the  state  space  approach  and  starting  with  the  displacement  equation 
(12).  They  have  conjectured  that  their  procedure  gives  rise  to  the  same  lossless  cascade 
model  as  the  one  obtained  via  the  Schur  algorithm.  However,  this  conjecture  was  never 
established. 

The  approach  of  Genin  et  al.  (which  is  essentially  the  same  as  that  of  Roberts  and 
Mullis)  involves  three  steps: 

(i)  embedding  eq.  (12)  into  a  state-space  model  {A,B,C,.D}  that  satisfies  the  loss¬ 
lessness  constraint  (11). 

(ii)  transforming  this  state-space  model  into  an  equivalent  balanced  form. 

(iii)  factoring  the  balanced  system  matrix  into  a  product  of  elementary  (i.e.  ,  single¬ 
state)  balanced  subsystem  matrices,  which  corresponds  to  a  cascade  decomposition 
of  the  system  described  by  {A,  B,  C,  D}. 


We  have  recently  constructed  an  alternative  synthesis  procedure  that  completely  avoids 
the  first  two  steps.  Our  procedure  determines  the  same  cascade  realization  as  the  one 
traditionally  obtained  by  preliminary  balancing,  but  it  does  so  by  (implicitly)  factoring 
the  unbalanced  model  {A,  B,  C,  D }  into  a  product  of  unbalanced  single-state  subsystems. 
Furthermore,  only  the  matrices  {A,B}  are  required  to  carry  out  our  procedure,  which 
makes  it  possible  to  avoid  the  embedding  step  (i). 


The  only  computation  involved  in  our  new  procedure  is  a  matrix  recursion,  viz. 


0 

\  G;+i  / 


=  G;  +  (Ti  -  I)  G{ 


PiJPi 


Wi 


(13) 


where  G0  =  B ,  /3,  is  the  first  row  of  ,  Wt  is  an  arbitrary  ./-unitary  matrix  and  Tt  is 
determined  (in  a  simple  manner)  by  A  alone  (Fig.  6). 
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Figure  6.  A  Cascade  Interpretation  of  the  Computational  Procedure  (13). 

This  recursion  is  the  matrix  equivalent  of  the  transfer  domain  recursion  (5),  which  was 
originally  derived  without  any  reference  to  lossless  state-space  models.  In  fact,  choosing 
A  =  L(/)  where  f(z)  :=  £f_0  /,-z*  makes  (13)  the  exact  equivalent  of  (5),  in  the  sense 
that  the  two  are  related  via  a  suitable  transform:  when  Q  =  0  for  all  i  this  is  the 
conventional  Z-transform,  but  for  non-vanishing  the  appropriate  transform  is  given  by 
a  Newton-series  expansion  (seeLev-Ari  and  Kailath  (1990);  Ackner,  Lev-Ari  and  Kailath 
(1990b)). 

The  relation  between  recursive  layer  peeling  and  the  factorization  of  matrices  is  now 

well  understood  (see  Section  2,  as  well  as  Lev-Ari  and  Kailath  (1984),  Bruckstein  and 

Kailath  (1987));  it  holds  even  for  matrices  with  the  generalized  displacement  representa¬ 
tion  (3),  viz., 

d(Z,  Z)R  =  GJG*  (14) 

where  G  denotes  the  matrix  of  coefficients  of  the  power-series  expansion  of  G(z),  (see, 
e.g.,  Lev-Ari  and  Kailath  (1986)).  Furthermore,  it  follows  that  for  every  d(z,w)  of  the 
form  (4)  and  for  every  (finite)  matrix  R, 

In  {d{ Z,  Z)R}  =  In  {d{ Z,  Z)R“^ }  (15) 

where  the  reversed  matrix  R-^  is  obtained  by  transposing  R-1  with  respect  to  the  anti¬ 
diagonal,  namely 

R-S  :=  I  (R-1)T  I  (16a) 
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where  the  superscript  T  denotes  the  conventional  (non-Hermitian)  transpose,  and  I  is 
the  anti-diagonal  unity  matrix,  viz., 


I  := 


0 


0 


v 


1 


(166) 


) 


The  fundamental  result  (15)  implies  that  R  and  R  "  have  the  same  displacement  struc¬ 
ture,  and  therefore  that  there  exists  a  matrix  H  such  that 


•rifZ.ZjlH  =  HJH 


(17) 


It  does  not  tell  us,  however,  how  to  obtain  H  or  what  is  its  corresponding  cascade  model. 
It  has  been  shown  in  a  limited  context  (and  by  fairly  lengthy  arguments)  that  the  cascade 
model  for  R-^  is  obtained  by  reversing  the  order  of  the  layers  in  the  cascade  model  for 
R  (Lev-Ari  and  Kailath  (1984)). 

One  of  the  advantages  of  the  state-space  approach  is  that  it  provides  a  direct  con¬ 
nection  between  the  representations  of  R  and  R-^.  Indeed  it  follows  directly  from  (11) 
that 


A 

BV 

(r- 

o  ^ 

1  A 

B  ) 

'  R-1 

o  N 

u 

w 

l  0 

J ) 

D) 

l  0 

J) 

Therefore  we  can  obtain  an  explicit  characterization  of  the  state-space  model  for  R 
viz., 

/  A  B  W  R-^  0  \  /  A  B  \  (  R-il  0  \ 

(18a) 


'  A  B  N 

U  D) 


\  o  J  ) 


\  C  DJ 


where 


'  A  B  ^ 

f*  °)| 

f  A  B  N 

X 

(  I  0  N 

U  Dj 

1 0  /  j  ' 

U  d  j 

lo  1  > 

(186) 
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Notice  that  A  =  AH  and  D  =  DT.  The  corresponding  transfer  function  is 

D  +  C(z~lI  -  A)“’B  =  Dt  +  B -  IAtI)-1ICt  =  [D  +  C (z^I  -  A)-1B]T 
and,  therefore, 

D  +  C(*-*I  -  A)-'B  =  0j(z) . . .  el(z)  (19) 

which  establishes  in  a  direct  manner  the  fact  that  the  cascade  model  for  R-^  is  obtained 
by  reversing  the  order  of  the  layers  (and  taking  the  transpose  of  each  ©1(2))  in  the  cascade 
model  for  R.  This  makes  it  possible  to  reconstruct  C  from  the  cascade  model  obtained 
for  {A,B}  and,  consequently,  to  obtain  explicit  expressions  for  R-1  ( i.e .  ,  Gohberg- 
Semencul  type  formulas)  as  well  as  solutions  of  linear  equations  involving  R.  In  fact, 
our  nev.  cascade  synthesis  method  gives  rise  to  two  new  algorithms  for  the  calculation 
of  C  for  all  matrices  with  a  displacement  structure.  In  the  particular  case  of  Toeplitz 
matrices  one  of  these  coincides  with  the  well-known  Levinson  algorithm,  while  the  other 
is  entirely  new. 

The  state-space  approach  has  also  been  extensively  used  in  the  solution  of  the  Nevanlinna- 
Pick  interpolation  problem,  which  arises,  among  other  applications  in  the  context  of 
model-order  reduction.  The  objective  of  model-order  reduction  is  to  approximate  a  given 
impulse  response  {hk}  by  a  rational  stable  transfer  function3 

=  (20) 

of  a  prescribed  order.  It  has  been  shown  that  when  the  quality  of  approximation  is 
measured  by  the  so-called  Hankel  norm  (applied  to  the  difference  between  hk  and  h^) 
then  hk  -  hk  is  the  causal  part  of  the  impulse  response  of  a  certain  non-causal  filter  that 
arises  in  the  following  interpolation  problem  (Genin  and  Rung  (1981)): 

Given  a  stable  and  causal  impulse  response  {hk}  find  a  function  f(z)  with  r 

poles  (or  less)  within  the  unit  disc  and  with  the  least  possible  infinity-norm 

3IIere  we  shall  maintain  the  convention  that  ©»(;)  is  a  first-order  matrix  polynomial  in  z,  while  H{z) 
is  nevertheless  defined  as  a  power  series  in  z-1. 
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|/|  |oo  such  that,  for  i  =  1, 2, . . . ,  n 


m  =  lira  {B(z)H(z)} 


where 

OO 

H(z)  =  £ 

fc=0 

is  the  transfer  function  associated  with  the  given  impulse  response,  {(,}  are 
the  poles  of  this  transfer  function  (all  of  which  are  within  the  unit  disc)  and 
B(z)  is  the  Blaschke  product  determined  by  these  poles,  i.e., 


b(z ) = n 

3 


1-Qz 


The  desired  reduced-order  transfer  function  H(z)  is  obtained  via 


(21) 


H(z)  =  H(z)  —  causal  part  of 


(22) 


If  H(z)  is  a  matrix  function  (corresponding  to  a  MIMO  system)  then  so  are  H(z)  and 
f(z );  in  this  case  {£}  are  the  Smith-MacMillan  poles  of  H(z). 


Though  f(z)  can  be  determined  via  the  Nevanlinna-Pick  procedure,  it  has  become 
customary  to  solve  this  problem  via  the  state  space  approach.  Starting  with  a  state- 
space  model  {A,B,C,D}  for  the  given  impulse  response  one  obtains  first  a  balanced 
realization  for  the  same  transfer  function  and  then,  applying  a  transformation  described 
by  Glover  (1984),  determines  the  state-space  model  {A,B,C,  D}  of  the  reduced-order 
transfer  function  H(z). 


The  preference  for  working  in  the  state-space  domain  has  been,  at  least  partly,  mo¬ 
tivated  by  the  lack  of  efficient  computational  procedures  fur  the  matrix  Nevanlinna-Pick 
problem.  The  standard  procedure  for  solving  this  problem  (see,  e.g.,  Delsarte,  Genin, 
and  Kamp  (1979))  requires  computationally-intensive  matrix  operations.  The  tangential 
Schur  algorithm  of  Dewilde  and  Dym  allows  a  significant  reduction  in  computational 
requirements  (though,  strictly  speaking  it  does  not  apply  to  the  case  of  functions  with 
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poles  within  the  unit  circle).  A  further  simplification  would  be  achieved  if  our  tangential 
Schur  algorithm  can  be  generalized  to  arbitrary  {£,};  currently  it  applies  only  to  £»  =  0 
for  all  i.  Since  our  procedure  involves  only  scalar  operations  and  determines  the  matrix 
function  f(z)  directly  from  the  given  impulse  response  H(z),  it  can  provide  an  attrac¬ 
tive  alternative  to  the  current  state-space  method  for  solving  the  model-order  reduction 
problem. 

5.  WAVELET  REPRESENTATION  OF  RANDOM  PROCESSES 

Wavelet  representations  of  deterministic  functions  have  recently  become  an  area  of  active 
research  in  a  variety  of  disciplines.  Applications  have  been  found  in  the  fields  of  Math¬ 
ematics  (Mallat  (1989)),  Physics  (Daubechies,  Grossman  and  Meyer  (1986)),  Numerical 
Analysis  (Glovvinski  et  al.  (1990))  and  Signal  Processing  (Mallat  (1989)).  However,  rel¬ 
atively  few  results  are  available  for  the  properties  of  wavelet  representations  of  random 
processes,  the  primary  difficulty  being  that  the  existing  theory  of  wavelets  has  been  ex¬ 
clusively  developed  in  the  context  of  Z2(IR).  We  have  shown  that  wavelet  representations 
can  be  extended  to  random  processes  with  either  finite  energy  or  finite  power.  We  have 
also  shown  that  periodically-correlated  processes,  namely  processes  with  periodic  statis¬ 
tics,  give  rise  to  wavelet  representations  with  periodic  correlation  structures  (Genossar, 
Goldburg,  Lev-Ari  and  Kailath  (1990)). 

The  wavelet  representation  at  resolution  level  l  of  a  function  /  €  T2(1R)  is 

fl(t)  =  ao,kKk{t)  iEE  6m,fctfm1fc(0>  1  >  0  >  (23a) 

A-eZ  m=0AgZ 

where  ip(t)  is  a  wavelet  function  and  <f>(t)  is  a  closely  related  function  known  as  the  scaling 
function  associated  with  the  wavelet  0(i).  Dilations  and  shifts  of  these  two  functions, 
viz., 

Mi)--=2',2Wt-kP)  ,  Mt)  :=  2"V(2'(  -  kP),  (236) 

give  rise  to  a  family  of  mutually  orthonormal  functions.  Consequently,  the  coefficients  of 
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the  representation  (23a)  are  given  by  inner  products,  viz., 


k,k  ■=  and  a0ik:=  (f,<j>0,k)- 


(23  c) 


Convergence  results 

We  have  shown  (in  Genossar,  Goldburg,  Lev-Ari,  and  Kailath  (1990))  that  for  wavelets 
with  finite  support  the  representation  at  resolution  level  l  is  well  defined  even  for  certain 
processes  that  are  not  members  of  Z2(1R)  and,  in  particular,  for  processes  with  finite 
power.  Moreover,  the  /- th  level  approximation  ffit)  converges  to  the  original  function  in 
the  following  sense:  For  any  compact  interval  I  C  IR, 

J  —  — *  0»  as  /  —»  oo  .  (24) 

For  wavelet  representations  of  random  processes ,  whether  of  finite  energy  or  of  finite 
power,  we  use  the  same  definition  as  equation  (23).  Here  too,  we  have  shown  that  the 
coefficients  and  approximations  are  well  defined.  The  corresponding  convergence  results 
are  expressed  in  terms  of  mean  square  errors.  For  random  processes  with  finite  energy 
we  have  shown  that 

Hm  J  £  ||x(t)  —  .i{(t)|2|  dt  =  0  .  (25) 

For  wavelet  representations  of  random  processes  with  finite  power  (and  using  wavelets  of 
compact  support)  we  have  shown  that  for  any  compact  interval  / 

Hm  £  ||a:(t)  —  £j(t)|2}  dt  =  0  .  (26) 

Periodically-Correlated  Processes 

Periodically-correlated  processes  provide  a  particularly  interesting  example  of  wavelet 
representations  for  random  processes.  A  periodically-correlated  process  with  period  P 
(also  referred  to  as  a  cyclostationary  process)  is  a  random  process  {x(t) ;  — oo  <  t  <  oo} 


22 


with  finite  second  moments,  whose  mean  and  autocovariance  functions  satisfy  the  condi¬ 
tions: 

£{x{t)}  =  £{x(t  +  P)}  (27  a) 

r(t,  s)  :=  £  (a:(i)a:*(s)}  =  r(t  -f  P,  $  +  P)  (276) 

The  wavelet  representation  of  this  process  gives  rise  to  a  filter-bank  model  (Fig.  7). 


Figure  7.  A  Filter  Model  for  Periodically-Correlated  Processes 

The  summation  of  the  outputs  of  the  2l  wavelet  filters  at  resolution  level  /  yields  xj(t), 
the  detail  signal  at  level  l.  We  have  shown  that  if  x(t)  is  a  periodically-correlated  process 
with  period  P,  then  all  input  sequences  to  this  filter-bank  model,  namely  the  sequence 
«o  ,k  as  well  as  the  sequences  6<t2n+r,  0</,  0<r<2i— 1  (for  fixed  /,  r),  are  jointly  wide- 
sense-stationary  in  the  index  k  (see  Genossar,  Goldburg,  Lev-Ari,  and  Kailath  (1990)). 
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